# =====================================================================
# Population Density Declines by Distance to a Business District
# ---------------------------------------------------------------------
# Input : Processed_Data/Intermediate/Section_2/distance_matrix/<City>_<comb>.rds
# Output: Processed_BD_Analysis/<City>_<comb>.rds
#         Results/Table_1/Table_1_<comb>.xlsx
#         Results/Table_1/Table_A1_<comb>.xlsx
# =====================================================================

rm(list = ls())
suppressPackageStartupMessages({
  library(data.table)
  library(fixest)
  library(modelsummary)
  library(tibble)
  library(clubSandwich)
  library(here)
})

# --- Paths (here-based) ------------------------------------------------
ROOT    <- normalizePath(here::here(), winslash = "/")
in_dir  <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2", "distance_matrix")
out_dir <- file.path(ROOT, "Processed_Data", "Intermediate", "Section_2", "Processed_BD_Analysis")
res_dir <- file.path(ROOT, "Results", "Table_1")

dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(res_dir, recursive = TRUE, showWarnings = FALSE)

# --- Parameters -------------------------------------------------------
combinations <- c(
  "80K_04_2020","80K_09_2020","80K_14_2020",
  "90K_04_2020","90K_09_2020","90K_14_2020",
  "100K_04_2020","100K_09_2020","100K_14_2020",
  "110K_04_2020","110K_09_2020","110K_14_2020"
)

cities <- c("Johannesburg","Capetown","Tshwane","Durban","NM_Bay",
            "Nairobi","Lagos","Abidjan","Cairo","London","Istanbul","Seoul",
            "SaoPaulo","KualaLumpur","Moscow","Mexico_City","Buenos_Aires",
            "Shenzhen","Frankfurt","Zurich","Bogota","Tehran","Dubai",
            "New_York","Mumbai")

# --- Skip logic -------------------------------------------------------
should_skip <- function(city, comb) {
  C <- paste0("Comb_", comb)
  if ((city %in% c("Lagos","Mumbai") && C %in% c("Comb_100K_14_2020","Comb_110K_14_2020")) ||
      (city == "Mumbai" && C == "Comb_90K_14_2020") ||
      (city %in% c("Lagos","Buenos_Aires","Mumbai") && C == "Comb_110K_09_2020") ||
      (city %in% c("Abidjan","Buenos_Aires") && C == "Comb_110K_14_2020")) TRUE else FALSE
}

# =====================================================================
# MAIN LOOP
# =====================================================================
for (comb in combinations) {
  
  # -------------------------------------------------------------------
  # Phase 1: process each city
  # -------------------------------------------------------------------
  for (city in cities) {
    
    if (should_skip(city, comb)) next
    
    infile  <- file.path(in_dir,  paste0(city, "_", comb, ".rds"))
    outfile <- file.path(out_dir, paste0(city, "_", comb, ".rds"))
    if (!file.exists(infile)) { message("Missing: ", infile); next }
    
    dt <- readRDS(infile) |> as.data.table()
    
    # Expected columns: x, y, building, pop, ntl, min_dist_cbd, bd_dummy, building_cutoff, cluster_size, year, city
    dt <- dt[!is.na(ntl) & !is.na(pop) & pop >= 10]
    
    # Distance bins
    dt[, bounds := fcase(
      min_dist_cbd <  5000,  5,
      min_dist_cbd < 10000, 10,
      min_dist_cbd < 15000, 15,
      min_dist_cbd < 20000, 20,
      min_dist_cbd < 25000, 25,
      default = 30
    )]
    
    dt[, counts := .N, by = bounds]
    
    # Collapse to bounds
    agg <- dt[, .(
      pop_sum = sum(pop, na.rm = TRUE),
      ntl_sum = sum(ntl, na.rm = TRUE),
      counts  = max(counts, na.rm = TRUE)
    ), by = bounds]
    
    agg[, pop_density := pop_sum / counts]
    agg[, ntl_density := ntl_sum / counts]
    agg[, pop_km_density := NA_real_]
    
    # Reference within 5 km
    pop5  <- agg[bounds == 5, pop_sum]
    popd5 <- agg[bounds == 5, pop_density]
    ntld5 <- agg[bounds == 5, ntl_density]
    
    agg[, population_5    := ifelse(length(pop5),  max(pop5,  na.rm = TRUE), NA_real_)]
    agg[, rel_pop_5       := pop_sum / population_5]
    agg[, population_d_5  := ifelse(length(popd5), max(popd5, na.rm = TRUE), NA_real_)]
    agg[, rel_pop_d_5     := pop_density / population_d_5]
    agg[, rel_pop_diff_5  := (pop_density - population_d_5) * 100 / population_d_5]
    agg[, ntl_d_5         := ifelse(length(ntld5), max(ntld5, na.rm = TRUE), NA_real_)]
    agg[, rel_ntl_d_5     := ntl_density / ntl_d_5]
    agg[, rel_ntl_diff_5  := (ntl_density - ntl_d_5) * 100 / ntl_d_5]
    agg[, city_name := city]
    
    saveRDS(agg, outfile)
  }
  
  # -------------------------------------------------------------------
  # Phase 2: combine cities
  # -------------------------------------------------------------------
  files_for_c <- file.path(out_dir, paste0(cities, "_", comb, ".rds"))
  files_for_c <- files_for_c[file.exists(files_for_c)]
  if (!length(files_for_c)) { message("No processed for ", comb); next }
  
  all_dt <- rbindlist(lapply(files_for_c, readRDS), use.names = TRUE, fill = TRUE)
  
  for (b in seq(5, 25, 5)) {
    all_dt[bounds == b, paste0("pop_b", b) := pop_sum]
    all_dt[, paste0("population_b", b) := max(get(paste0("pop_b", b))), by = city_name]
  }
  
  all_dt[, south_africa := city_name %in% c("Capetown","Johannesburg","Durban","Tshwane","NM_Bay")]
  all_dt[, other_africa := city_name %in% c("Lagos","Nairobi","Abidjan","Cairo")]
  
  all_dt[is.na(population_b20), population_b20 := 0]
  pop_by_city <- unique(all_dt[, .(city_name, population_b5, population_b10, population_b15, population_b20)])
  pop_by_city[, ratio_10_20 := (fifelse(is.na(population_b15), 0, population_b15) +
                                  fifelse(is.na(population_b20), 0, population_b20)) /
                (fifelse(is.na(population_b5),  0, population_b5) +
                   fifelse(is.na(population_b10), 0, population_b10))]
  all_dt <- merge(all_dt, pop_by_city[, .(city_name, ratio_10_20)], by = "city_name", all.x = TRUE)
  
  rename_map <- c("Capetown"="Cape Town","Mexico_City"="Mexico City",
                  "SaoPaulo"="São Paulo","KualaLumpur"="Kuala Lumpur",
                  "Tshwane"="Tshwane (Pretoria)","Durban"="eThekwini (Durban)",
                  "NM_Bay"="NM Bay (Port Elizabeth)",
                  "Buenos_Aires"="Buenos Aires","New_York"="New York City")
  all_dt[, city_name := ifelse(city_name %in% names(rename_map), rename_map[city_name], city_name)]
  
  glob <- all_dt[!south_africa & !other_africa,
                 .(global_rel_ntl_diff_5 = mean(rel_ntl_diff_5, na.rm = TRUE),
                   global_rel_pop_diff_5 = mean(rel_pop_diff_5, na.rm = TRUE)),
                 by = bounds]
  all_dt <- merge(all_dt, glob, by = "bounds", all.x = TRUE)
  
  all_dt[, diff_global_ntl := rel_ntl_diff_5 - global_rel_ntl_diff_5]
  all_dt[, diff_global_pop := rel_pop_diff_5 - global_rel_pop_diff_5]
  all_dt[, pop_km_density := pop_density * 100]
  all_dt[, ln_pop_km_density := log(pmax(pop_km_density, .Machine$double.eps))]
  all_dt[, city_code := as.integer(factor(city_name))]
  all_dt[, bounds_f := relevel(factor(bounds), ref = "20")]
  
  saveRDS(all_dt, file.path(out_dir, paste0("ALL_CITIES_", comb, ".rds")))
  
  # -------------------------------------------------------------------
  # Phase 3: regressions + formatted output
  # -------------------------------------------------------------------
  reg_dt <- all_dt[bounds < 25]
  reg_dt[, bounds_f := droplevels(bounds_f)]  # drop 25/30 levels
  
  # Models
  m1 <- feols(pop_km_density ~ bounds_f,                               data = reg_dt, cluster = ~city_code)
  m2 <- feols(pop_km_density ~ bounds_f * south_africa,                data = reg_dt, cluster = ~city_code)
  m3 <- feols(pop_km_density ~ bounds_f * south_africa | city_name,    data = reg_dt, cluster = ~city_code)
  m4 <- feols(pop_km_density ~ ntl_density + bounds_f * south_africa | city_name, data = reg_dt, cluster = ~city_code)
  
  m5 <- feols(ln_pop_km_density ~ bounds_f,                            data = reg_dt, cluster = ~city_code)
  m6 <- feols(ln_pop_km_density ~ bounds_f * south_africa,             data = reg_dt, cluster = ~city_code)
  m7 <- feols(ln_pop_km_density ~ bounds_f * south_africa | city_name, data = reg_dt, cluster = ~city_code)
  m8 <- feols(ln_pop_km_density ~ ntl_density + bounds_f * south_africa | city_name, data = reg_dt, cluster = ~city_code)
  
  # Labels
  coef_map <- c(
    "(Intercept)"                 = "Constant",
    "bounds_f5"                   = "0 km. < Distance to BD <= 5 km.",
    "bounds_f10"                  = "5 km. < Distance to BD < 10 km.",
    "bounds_f15"                  = "10 km. < Distance to BD < 15 km.",
    "south_africaTRUE"            = "South Africa",
    "bounds_f5:south_africaTRUE"  = "0–5 km. × S. Africa",
    "bounds_f10:south_africaTRUE" = "5–10 km. × S. Africa",
    "bounds_f15:south_africaTRUE" = "10–15 km. × S. Africa",
    "ntl_density"                 = "Avg. NTL Luminosity"
  )
  
  gof_map <- data.frame(
    raw   = c("nobs", "r.squared"),
    clean = c("Observations", "R-squared"),
    fmt   = c(0, 2)
  )
  
  country_dummy_row <- tibble(
    term  = "Country Dummy",
    `(1)` = "No", `(2)` = "No", `(3)` = "Yes", `(4)` = "Yes"
  )
  
  notes_row <- tibble(
    term  = "Notes",
    `(1)` = "*** p<0.01, ** p<0.05, * p<0.1. Robust t-statistics in brackets. SEs clustered by city. Base category: 15–20 km.",
    `(2)` = "", `(3)` = "", `(4)` = ""
  )
  
  # one file per combination in Results/Table_1/
  tab1_path <- file.path(res_dir, paste0("Table_1_",  comb, ".xlsx"))
  tabA1_path <- file.path(res_dir, paste0("Table_A1_", comb, ".xlsx"))
  
  # Table 1
  modelsummary(
    list(`(1)` = m1, `(2)` = m2, `(3)` = m3, `(4)` = m4),
    coef_map           = coef_map,
    statistic          = "([{statistic}])",
    statistic_vertical = TRUE,
    stars              = c('*' = .1, '**' = .05, '***' = .01),
    gof_map            = gof_map,
    add_rows           = rbind(country_dummy_row, notes_row),
    title              = "Table 1: Population Density Declines by Distance to a Business District",
    output             = tab1_path
  )
  
  # Appendix Table A1 (your log-density spec bundle)
  modelsummary(
    list(`(1)` = m5, `(2)` = m6, `(3)` = m7, `(4)` = m8),
    coef_map           = coef_map,
    statistic          = "([{statistic}])",
    statistic_vertical = TRUE,
    stars              = c('*' = .1, '**' = .05, '***' = .01),
    gof_map            = gof_map,
    add_rows           = rbind(country_dummy_row, notes_row),
    title              = "Table A1: Log Population Density Declines by Distance to a Business District",
    output             = tabA1_path
  )
  
  message("Completed ", comb)
}

message("All done.")


